I Предварительный анализ данных

if (interactive() && !str_ends(getwd(), "R/Statistical_Analysis_2022/CITY_US/Tkachenko")) {
    setwd("R/Statistical_Analysis_2022/CITY_US/Tkachenko")
}

data <- read_excel("../CITY_US.STD/CITY_shortname.xls")
data[data == "NA"] <- NA
data[, -(1:2)] <- data.frame(lapply(data[, -(1:2)], as.numeric))
fullnames <- names(read_excel("../CITY_US.STD/CITY.xls"))

1 Разобраться в том, что означают признаки.

print(fullnames)
##  [1] "CITY City"                                                                            
##  [2] "STATE State"                                                                          
##  [3] "AREA Land area, 1990 (sq.miles)"                                                      
##  [4] "R_AREA Land area, 1990 (sq.miles), ranks"                                             
##  [5] "POP92 Population,1992"                                                                
##  [6] "R_POP92 Population,1992, ranks"                                                       
##  [7] "POP80 Population,1980"                                                                
##  [8] "R_POP80 Population,1980, ranks"                                                       
##  [9] "POP_CH Population growth rate,1980-1992"                                              
## [10] "R_POP_CH Population growth rate,1980-1992, ranks"                                     
## [11] "POPDEN Population per square mile, 1992"                                              
## [12] "R_POPDEN Population per square mile, 1992, ranks"                                     
## [13] "OLD % older 65 years"                                                                 
## [14] "R_OLD % older 65 years, ranks"                                                        
## [15] "BLACK Black population,1990"                                                          
## [16] "R_BLACK Black population,1990, ranks"                                                 
## [17] "BLACK% % Black population,1990"                                                       
## [18] "R_BLACK% % Black population,1990, ranks"                                              
## [19] "ASIAN Asian or Pacific Islander population,1990"                                      
## [20] "R_ASIAN Asian or Pacific Islander population,1990, ranks"                             
## [21] "ASIAN% % Asian or Pacific Islander population,1990"                                   
## [22] "R_ASIAN% % Asian or Pacific Islander population,1990, ranks"                          
## [23] "HISP Hispanic population, 1990"                                                       
## [24] "R_HISP Hispanic population, 1990, ranks"                                              
## [25] "HISP% % Hispanic population, 1990"                                                    
## [26] "R_HISP% % Hispanic population, 1990, ranks"                                           
## [27] "BORN_F % persons of foreign born"                                                     
## [28] "R_BORN_F % persons of foreign born, ranks"                                            
## [29] "LANG % of persons, speaking language other than English at home, 1990"                
## [30] "R_LANG % of persons, speaking language other than English at home, 1990, ranks"       
## [31] "HH1 % 1-person households, 1990"                                                      
## [32] "R_HH1 % 1-person households, 1990, ranks"                                             
## [33] "FAMIL1 % one-parent family households, 1990"                                          
## [34] "R_FAMIL1 % one-parent family households, 1990, ranks"                                 
## [35] "DEATH Infant death rate per 1000 live births,1988"                                    
## [36] "R_DEATH Infant death rate per 1000 live births,1988, ranks"                           
## [37] "CRIME Crime rate per 100000 population, 1991"                                         
## [38] "R_CRIME Crime rate per 100000 population, 1991, rate"                                 
## [39] "SCHOOL % of elementary and high school enrollment in public schools, 1990"            
## [40] "R_SCHOOL % of elementary and high school enrollment in public schools, 1990, ranks"   
## [41] "DEGREE % of adults with a bachelor's degree or higher, 1990"                          
## [42] "R_DEGREE % of adults with a bachelor's degree or higher, 1990, ranks"                 
## [43] "INCOME Median household income, 1989"                                                 
## [44] "R_INCOME Median household income, 1989, ranks"                                        
## [45] "ASSIST % of households, receiving public assistance,1989"                             
## [46] "R_ASSIST % of households, receiving public assistance,1989, ranks"                    
## [47] "POVERT % of persons below poverty level, 1989"                                        
## [48] "R_POVERT % of persons below poverty level, 1989, ranks"                               
## [49] "OLB_BIL % of housing units built 1939 or earlier, 1990"                               
## [50] "R_OLD_BI % of housing units built 1939 or earlier, 1990, ranks"                       
## [51] "OWNER Median value of specified owner-occupied housing units ($), 1990"               
## [52] "R_OWNER Median value of specified owner-occupied housing units ($), 1990, ranks"      
## [53] "RENTER % renter-occupied housing units,1990"                                          
## [54] "R_RENTER % renter-occupied housing units,1990, ranks"                                 
## [55] "GROSS Median gross rent of specified renter-occupied housing units ($), 1990"         
## [56] "R_GROSS Median gross rent of specified renter-occupied housing units ($), 1990, ranks"
## [57] "CONDOM % condominimum occupied housing units, 1990"                                   
## [58] "R_CONDOM % condominimum occupied housing units, 1990, ranks"                          
## [59] "TRANSP % of workers using public transportation"                                      
## [60] "R_TRANSP % of workers using public transportation, ranks"                             
## [61] "UNEMP Unemployment rate, 1991"                                                        
## [62] "R_UNEMP Unemployment rate, 1991, ranks"                                               
## [63] "LABOR Labor force - % change 1980-1990"                                               
## [64] "R_LABOR Labor force - % change 1980-1990, ranks"                                      
## [65] "LAB_F Female civilian labor force participation rate, 1990"                           
## [66] "R_LAB_F Female civilian labor force participation rate, 1990, ranks"                  
## [67] "MANLAB % of civilian labor force emploied in manufacturing, 1990"                     
## [68] "R_MANLAB % of civilian labor force emploied in manufacturing, 1990, ranks"            
## [69] "TAXE City government taxes per capita, 1991"                                          
## [70] "R_TAXE City government taxes per capita, 1991, ranks"                                 
## [71] "TEMPER Average daily temperature in July"                                             
## [72] "R_TEMPER Average daily temperature in July, ranks"                                    
## [73] "PRECEP Annual precipitation"                                                          
## [74] "R_PRECEP Annual precipitation"

2 Отобрать признаки

names_interesting <- c("AREA", "POP80", "POP92", "POPDEN", "CRIME", "BORN_F", "POVERT", "INCOME", "UNEMP", "TEMPER")

data <- data %>% select(all_of(c("CITY", "STATE", names_interesting)))

print(head(data))
## # A tibble: 6 × 12
##   CITY  STATE  AREA  POP80  POP92 POPDEN CRIME BORN_F POVERT INCOME UNEMP TEMPER
##   <chr> <chr> <dbl>  <dbl>  <dbl>  <dbl> <dbl>  <dbl>  <dbl>  <dbl> <dbl>  <dbl>
## 1 NEW_… NY     309. 7.07e6 7.31e6  23671  9236   28.4   19.3  29823   8.6   76.8
## 2 LOS_… CA     469. 2.97e6 3.49e6   7436  9730   38.4   18.9  30925   9     74.3
## 3 CHIC… IL     227. 3.01e6 2.77e6  12185    NA   16.9   21.6  26301   8.4   75.1
## 4 HOUS… TX     540. 1.60e6 1.69e6   3131 10824   17.8   20.7  26261   6.1   83.5
## 5 PHIL… PA     135. 1.69e6 1.55e6  11492  6835    6.6   20.3  24603   8     76.7
## 6 SAN_… CA     324  8.76e5 1.15e6   3546  8537   20.9   13.4     10   6.2   71

3 Определить вид признаков

Город и штат качественные, остальные количественные, ранги были порядковыми.

find_mode_freq <- function(x) {
    x <- x[!is.na(x)]
    return(max(tabulate(match(x, x))))
}

print(data %>% summarise(across(
    all_of(names_interesting),
    find_mode_freq
)))
## # A tibble: 1 × 10
##    AREA POP80 POP92 POPDEN CRIME BORN_F POVERT INCOME UNEMP TEMPER
##   <int> <int> <int>  <int> <int>  <int>  <int>  <int> <int>  <int>
## 1     1     1     1      2     1      4      2      1     5      3
print(sort(data$UNEMP))
##  [1]  2.3  3.2  3.8  3.9  4.1  4.4  4.6  4.7  4.7  4.8  4.8  5.0  5.0  5.0  5.0
## [16]  5.1  5.1  5.1  5.3  5.4  5.4  5.4  5.4  5.4  5.6  5.7  5.9  5.9  5.9  6.0
## [31]  6.0  6.1  6.1  6.1  6.1  6.2  6.2  6.4  6.4  6.5  6.7  6.7  6.7  6.8  6.9
## [46]  6.9  7.0  7.1  7.2  7.2  7.3  7.3  7.3  7.5  7.6  7.7  7.7  7.7  8.0  8.1
## [61]  8.4  8.4  8.5  8.6  9.0  9.0  9.4  9.5  9.6 10.0 10.4 10.6 10.7 11.0 11.9
## [76] 12.6 13.1

Все количественные буду считать непрерывными. возможно UNEMP непрерывный с плохой точностью.

4 не актуально

5 Построить matrix plot

if (interactive()) pdf("ggpairs_unedited.pdf")
ggpairs(
    data[, -(1:2)],
    lower = list(continuous = wrap("points", alpha = 0.5, size = 0.3)),
    diag = list(continuous = "barDiag")
)

if (interactive()) dev.off()

7 outliers

Убираю outliers:

Помечаю некорректные данные в INCOME как NA. Удаляю город из Аляски за плотность населения. Флорида выделсется на BORN_F-INCOME Гаваи выделяются низким уровнем безработицы. странно?

data$INCOME[data$INCOME < 100] <- NA
data <- data %>% filter(STATE != "AK")

6 Несимметричные распределения

Функция, которая логарифмирует, если это сделает выборку симметричнее

log_asymmetric <- function(x) {
    if (skewness(x, na.rm = TRUE) < abs(skewness(log(x), na.rm = TRUE))) {
        print("default")
        return(x)
    } else {
        print("logged")
        return(log(x))
    }
}

Автоматически логарифмирую то что имеет асимметрию и длинный хвост справа

data_logged <- data %>%
    mutate(across(all_of(names_interesting), log_asymmetric))
## [1] "logged"
## [1] "logged"
## [1] "logged"
## [1] "logged"
## [1] "logged"
## [1] "logged"
## [1] "default"
## [1] "default"
## [1] "logged"
## [1] "default"
if (interactive()) pdf("ggpairs_logged.pdf")
ggpairs(
    data_logged[, -(1:2)],
    lower = list(continuous = wrap("points", alpha = 0.5, size = 0.3)),
    diag = list(continuous = "barDiag")
)

if (interactive()) dev.off()

8 однородность

выглядит однородно.

9 не актуально

10 всякие характеристики

print_characteristics <- function(x) {
    list(
        mean = mean(x, na.rm = TRUE),
        var = var(x, na.rm = TRUE),
        skewness = skewness(x, na.rm = TRUE),
        kurtosis = kurtosis(x, na.rm = TRUE) - 3
    )
}

sapply(data_logged %>% select(all_of
(names_interesting)), print_characteristics)
##          AREA       POP80     POP92     POPDEN     CRIME      BORN_F    
## mean     4.754638   12.9069   13.01226  8.257586   9.206558   1.959747  
## var      0.6658984  0.5142076 0.4464288 0.5099096  0.06798296 0.8435754 
## skewness 0.1317646  1.453879  1.743773  0.1015321  0.1476689  0.3081197 
## kurtosis -0.3907363 3.033026  3.885752  -0.1531417 0.08872069 -0.7431516
##          POVERT     INCOME     UNEMP      TEMPER    
## mean     18.11842   25282.48   1.873674   77.60132  
## var      34.68872   14321837   0.09967048 37.59853  
## skewness 0.2243883  -0.2196584 -0.2186434 -0.1426746
## kurtosis -0.3539569 -0.6467894 0.6594639  0.7472063

II О виде распределений и о сравнении распределений

1 выполнить первое второе задание в логичном порядке

2 анализ вида распределения признаков

2. Описать распределения признаков визуально (с помощью нормальной бумаги и PP-plot) и по критерию хи-квадрат и другим критериям.

Сюда входит: normal probability plot (что это такое?), проверка по критериям Лиллиефорса, AD, хи-квадрат, Шапиро-Уилка. По критерию хи-квадрат, а также визуально по PP-plot можно проверить и гипотезы о согласии с другими распределениями, например, логнормальным.

Рассматриваем нормальность

все плотности.

if (interactive()) pdf(file = "all_density_plots.pdf")
ggplot(melt(data_logged, id = c("CITY", "STATE")), aes(value)) +
    theme_bw() +
    geom_histogram(aes(y = ..density..), bins = 15) +
    # geom_density() +
    facet_wrap(~variable, scales = "free") +
    geom_line(aes(y = dnorm(value,
        mean = tapply(value, variable, mean, na.rm = TRUE)[PANEL],
        sd = tapply(value, variable, sd, na.rm = TRUE)[PANEL]
    )), color = "red") +
    theme(legend.position = "bottom") +
    labs(x = "", y = "")

if (interactive()) dev.off()

Normal probability plot

if (interactive()) pdf(file = "normal_probability_plot.pdf")
ggplot(melt(data_logged, id = c("CITY", "STATE")), aes(sample = value)) +
    stat_qq_point(size = 2) +
    geom_abline() +
    facet_wrap(~variable, scales = "free")

if (interactive()) dev.off()

PP-plot

if (interactive()) pdf(file = "PP_plot.pdf")
ggplot(melt(data_logged, id = c("CITY", "STATE")), aes(sample = value)) +
    stat_pp_point(size = 2) +
    facet_wrap(~variable, scales = "free") +
    stat_pp_line()

if (interactive()) dev.off()

Всякие тесты на нормальность

data_logged$F_TEMPER <- NULL
test_all <- function(data, test, ...) {
    print(test(1:10, ...)$method)
    print(data %>% summarise(across(
        everything(),
        function(x) {
            test(x, ...)$p.value
        }
    )))
}

test_all(data_logged[, -(1:2)], shapiro.test)
## [1] "Shapiro-Wilk normality test"
## # A tibble: 1 × 10
##    AREA     POP80       POP92 POPDEN CRIME BORN_F POVERT INCOME UNEMP TEMPER
##   <dbl>     <dbl>       <dbl>  <dbl> <dbl>  <dbl>  <dbl>  <dbl> <dbl>  <dbl>
## 1 0.757 0.0000113 0.000000310  0.673 0.725 0.0892  0.675  0.463 0.487  0.368
test_all(data_logged[, -(1:2)], lillie.test)
## [1] "Lilliefors (Kolmogorov-Smirnov) normality test"
## # A tibble: 1 × 10
##    AREA  POP80    POP92 POPDEN CRIME BORN_F POVERT INCOME UNEMP TEMPER
##   <dbl>  <dbl>    <dbl>  <dbl> <dbl>  <dbl>  <dbl>  <dbl> <dbl>  <dbl>
## 1 0.536 0.0118 0.000633  0.291 0.580  0.148  0.646  0.866 0.746  0.766
test_all(data_logged[, -(1:2)], ad.test)
## [1] "Anderson-Darling normality test"
## # A tibble: 1 × 10
##    AREA    POP80       POP92 POPDEN CRIME BORN_F POVERT INCOME UNEMP TEMPER
##   <dbl>    <dbl>       <dbl>  <dbl> <dbl>  <dbl>  <dbl>  <dbl> <dbl>  <dbl>
## 1 0.492 0.000167 0.000000752  0.370 0.620 0.0362  0.830  0.692 0.697  0.371
test_all(data_logged[, -(1:2)], pearson.test)
## [1] "Pearson chi-square normality test"
## # A tibble: 1 × 10
##    AREA  POP80  POP92 POPDEN CRIME BORN_F POVERT INCOME UNEMP TEMPER
##   <dbl>  <dbl>  <dbl>  <dbl> <dbl>  <dbl>  <dbl>  <dbl> <dbl>  <dbl>
## 1 0.697 0.0377 0.0179  0.664 0.728  0.265  0.729  0.857 0.897 0.0626

Критерии отвергают нормальность распределения населения. В остальных признаках отклонения от нормального (логнормального для ранее логарифмированных признаков) распределения тестами не обнаружено.

1. Описать разницу между численностью населения в 1980 и 1992 годах визуально (с помощью ящиков с усами) и по критериям. Добавить признак, который разбивает города на “холодные” и “жаркие”. Сравнить их по характеристикам.

добавим признак

temper_split <- median(data_logged$TEMPER)
data_logged$F_TEMPER <- factor(ifelse(data_logged$TEMPER < temper_split, "cold", "hot"))

3 Сначала имеет смысл посмотреть на сравнение распределений в группах с помощью ящиков с усами

if (interactive()) pdf("boxplot_pop.pdf")
ggplot(
    melt(data_logged, id.vars = "F_TEMPER", measure.vars = c("POP80", "POP92")),
    aes(x = variable, y = value, color = F_TEMPER)
) +
    geom_boxplot() +
    scale_color_manual(values = c("blue", "red"))

if (interactive()) dev.off()


if (interactive()) pdf("boxplot_pop_all.pdf")
ggplot(
    melt(data_logged, id.vars = "F_TEMPER", measure.vars = c("AREA", "POP80", "POPDEN", "CRIME", "BORN_F", "POVERT", "INCOME", "UNEMP", "TEMPER")),
    aes(x = variable, y = value, color = F_TEMPER)
) +
    geom_boxplot() +
    facet_wrap(~variable, scales = "free") +
    scale_color_manual(values = c("blue", "red"))
## Warning: Removed 11 rows containing non-finite values (`stat_boxplot()`).

if (interactive()) dev.off()

Сравнение признаков при разбиении городов на группы по температуре

4 t-test

Разделение на группы было сделано по медиане, поэтому дизайн сбалансированный, и для t-теста можно не выяснять равенство дисперсий.

Но из интереса применю тест Фишера, а перед этим проверю группы на нормальность распределения.

Признак CRIME

Из ящиков с усами ожидаю отвержение нулевых гипотез в тестах.

hot <- (data_logged %>%
    filter(F_TEMPER == "hot"))$CRIME
cold <- (data_logged %>%
    filter(F_TEMPER == "cold"))$CRIME
cat(shapiro.test(hot)$p.value, shapiro.test(cold)$p.value)
## 0.5012446 0.954743
cat(pearson.test(hot)$p.value, pearson.test(cold)$p.value)
## 0.1314607 0.4171623

Тесты на нормальность не обнаружили ненормальности в распределениях. Применю тест Фишера

var.test(hot, cold)
## 
##  F test to compare two variances
## 
## data:  hot and cold
## F = 1.6011, num df = 37, denom df = 36, p-value = 0.1608
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.8269838 3.0904136
## sample estimates:
## ratio of variances 
##           1.601064

Тест Фишера не отверг гипотезу о равенстве дисперсий, но это не так важно, потому что в t-test дизайн сбалансирован и формула не изменится.

t.test(cold, hot)
## 
##  Welch Two Sample t-test
## 
## data:  cold and hot
## t = -3.2207, df = 70.063, p-value = 0.00194
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.29506315 -0.06938105
## sample estimates:
## mean of x mean of y 
##  9.114233  9.296455

t-test уверенно отвергает гипотезу о равенстве матожиданий, как и ожидалось.

5 непараметрический тест Вилкоксона.

Если бы дизайн был несбалансирован или не было бы нормальности распределений и нас не удовлетворяла асимптотичность t-критерия или нас интересовала другая характеристика положение, то мы бы применили критерий Вилкоксона.

wilcox.test(cold, hot, paired = FALSE)
## 
##  Wilcoxon rank sum exact test
## 
## data:  cold and hot
## W = 413, p-value = 0.001867
## alternative hypothesis: true location shift is not equal to 0

Он тоже отвергает схожесть распределений.

7 критерий Колмогорова-Смирнова

ks.test(cold, hot)
## 
##  Exact two-sample Kolmogorov-Smirnov test
## 
## data:  cold and hot
## D = 0.38762, p-value = 0.004648
## alternative hypothesis: two-sided

Тест Колмогорова-Смирнова обнаруживает разницу в форме распределений.

Признак INCOME

для этого признака многое повторится, но я ожидаю неотвержения нулевых гипотез в тестах, потому что ящики с усами похожи.

Сейчас я буду применять t-test для сравнения признака INCOME в группах разделения по температуре. У меня сбалансированный дизайн (в группах по 38 городов), поэтому не важно, равны дисперсии или нет, но я проверю. Применять тест Фишера для проверки равенства дисперсий можно только для нормально распределенных (внутри групп) признаков.

hot <- (data_logged %>%
    filter(F_TEMPER == "hot"))$INCOME
cold <- (data_logged %>%
    filter(F_TEMPER == "cold"))$INCOME

Проверяю нормальность в группах тестом Пирсона (хи квадрат) и Шапира-Уилка

cat(shapiro.test(hot)$p.value, shapiro.test(cold)$p.value)
## 0.8588537 0.3617601
cat(pearson.test(hot)$p.value, pearson.test(cold)$p.value)
## 0.7321958 0.4158802

Тесты не обнаружили ненормальности в распределениях в группах. Поэтому применю тест Фишера.

var.test(hot, cold)
## 
##  F test to compare two variances
## 
## data:  hot and cold
## F = 1.0373, num df = 33, denom df = 31, p-value = 0.9211
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.5101231 2.0928364
## sample estimates:
## ratio of variances 
##           1.037312

Тест Фишера не отвергает гипотезу о равенстве дисперсий.

Наконец применю t-test, никак не используя результаты про наверно (гипотезу нельзя принимать) равенство дисперсий, потому что дизайн сбалансирован. При нормальных распределениях тест точный.

t.test(cold, hot)
## 
##  Welch Two Sample t-test
## 
## data:  cold and hot
## t = 0.16315, df = 63.88, p-value = 0.8709
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1722.032  2028.311
## sample estimates:
## mean of x mean of y 
##  25361.38  25208.24

t-test не отвергает гипотезу о равенстве матожиданий. Применю тест Вилкоксона. Он ранговый и проверят другую гипотезу про сумму рангов.

wilcox.test(cold, hot)
## 
##  Wilcoxon rank sum exact test
## 
## data:  cold and hot
## W = 569, p-value = 0.7548
## alternative hypothesis: true location shift is not equal to 0

Тест Вилкоксона не обнаружил разницы в распределениях, так как не отверг нулевую гипотезу.

В тесте Колмогорова-Смирнова есть смысл, если прошлые тесты не нашли различий. Тест Колмогорова-Смирнова сравнивает функции распределения целиком и проверяет, похожи ли формы распределений.

ks.test(cold, hot)
## 
##  Exact two-sample Kolmogorov-Smirnov test
## 
## data:  cold and hot
## D = 0.15993, p-value = 0.7266
## alternative hypothesis: two-sided

p-value большое, тест не нашел заметной разницы между формами распределений.

8 сравнение зависимых выборок

Применяю t-test для зависимых выборок.

t.test(
    data_logged$POP80,
    data_logged$POP92,
    paired = TRUE
)
## 
##  Paired t-test
## 
## data:  data_logged$POP80 and data_logged$POP92
## t = -4.7135, df = 75, p-value = 1.098e-05
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.14988791 -0.06083053
## sample estimates:
## mean difference 
##      -0.1053592

Благодаря тому что тест парный, он смог отвергнуть гипотезу о равенстве матожиданий.

Статистика зависимого критерия больше статистики независимого (при положительной корреляции), поэтому мощнее.

Парный тест Вилкоксона

wilcox.test(data_logged$POP80, data_logged$POP92, paired = TRUE)
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  data_logged$POP80 and data_logged$POP92
## V = 730, p-value = 0.0001492
## alternative hypothesis: true location shift is not equal to 0

Из-за того что тесты парные, они обнаруживают разницу в распределениях.

III Об анализе зависимостей

1 вспомним pairs plot

if (interactive()) pdf("ggpairs_logged.pdf")
ggpairs(
    data_logged[, -(1:2)],
    lower = list(continuous = wrap("points", alpha = 0.5, size = 0.3)),
    diag = list(continuous = "barDiag")
)

if (interactive()) dev.off()
if (interactive()) pdf("ggpairs_logged_colored.pdf")
ggpairs(
    data_logged[, -(1:2)],
    aes(color = F_TEMPER),
    lower = list(continuous = wrap("points", alpha = 0.5, size = 0.6)),
    diag = list(continuous = "barDiag"),
    columns = names_interesting
) +
    scale_color_manual(values = c("blue", "red")) +
    scale_fill_manual(values = c("blue", "red"))

if (interactive()) dev.off()

2 Начинать нужно с анализа линейных зависимостей

Матрица корреляций пирсона

data_logged <- data_logged %>%
    relocate(TEMPER, AREA, POP92, POP80, POPDEN, BORN_F, UNEMP, POVERT, CRIME, INCOME)


if (interactive()) pdf(file = "cor_matrix_pearson.pdf")
ggplot(melt(cor(data_logged %>% select(-CITY, -STATE, -F_TEMPER), method = "pearson", use = "pairwise.complete.obs"))) +
    geom_raster(aes(x = Var2, y = Var1, fill = value)) +
    geom_text(aes(x = Var2, y = Var1, label = round(value, 2))) +
    scale_fill_gradient2() +
    theme_dark() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

if (interactive()) dev.off()

3 Затем можно переходить к ранговым коэффициентам корреляции.

Матрица корреляций спирмана

if (interactive()) pdf(file = "cor_matrix_spearman.pdf")
ggplot(melt(cor(data_logged %>% select(-CITY, -STATE, -F_TEMPER), method = "spearman", use = "pairwise.complete.obs"))) +
    geom_raster(aes(x = Var2, y = Var1, fill = value)) +
    geom_text(aes(x = Var2, y = Var1, label = round(value, 2))) +
    scale_fill_gradient2() +
    theme_dark() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

if (interactive()) dev.off()

4 Проинтерпретируйте найденные корреляции - можно ли сказать, что является причиной, что следствием

частные корреляции

Вычитаю из иностранных и плотности доход

(data_logged %>%
    dplyr::select(BORN_F, POPDEN, INCOME) %>%
    drop_na() %>%
    pcor(method = "pearson"))$estimate["BORN_F", "POPDEN"]
## [1] 0.5507378
(data_logged %>%
    dplyr::select(BORN_F, POPDEN, TEMPER) %>%
    drop_na() %>%
    pcor(method = "pearson"))$estimate["BORN_F", "POPDEN"]
## [1] 0.5038904

Частная корреляция иностранно рожденных и плотности населения за вычетом дохода больше чем обычная корреляция.

if (interactive()) pdf(file = "pcor.pdf")

ggplot(melt((data_logged %>%
    dplyr::select(UNEMP, POVERT, INCOME, CRIME) %>%
    drop_na() %>%
    pcor(method = "pearson"))$estimate), aes(x = Var2, y = Var1)) +
    geom_raster(aes(fill = value)) +
    geom_text(aes(label = round(value, 2))) +
    scale_fill_gradient2()

if (interactive()) dev.off()

if (interactive()) pdf(file = "cor.pdf")
ggplot(melt(cor(data_logged %>%
    dplyr::select(UNEMP, POVERT, INCOME, CRIME) %>%
    drop_na())), aes(x = Var2, y = Var1)) +
    geom_raster(aes(fill = value)) +
    geom_text(aes(label = round(value, 2))) +
    scale_fill_gradient2()

if (interactive()) dev.off()


if (interactive()) pdf(file = "pairs.pdf")

ggpairs(
    data_logged %>%
        dplyr::select(UNEMP, POVERT, INCOME, CRIME) %>%
        drop_na(),
    lower = list(continuous = wrap("points", alpha = 0.5, size = 0.3)),
    diag = list(continuous = "barDiag")
)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

if (interactive()) dev.off()